%% bdary_node (for Dirichlet)
function bdn = bdary_node(mesh, fespace, index)
bdf = bdary_face(mesh, index);
bdn = unique(bdf(:,2:4), "stable");
if "P2" == fespace.basis_type
    tri = bdf(:,2:4);
    t = [1 2;
        2 3;
        3 1];
    edge = reshape(tri(:,t), [], 2);
    edge = sort(edge,2);
    idx = zeros(size(edge,1),1);
    for e = 1:size(edge,1)
        key = sprintf("%d %d", edge(e,1), edge(e,2));
        idx(e) = mesh.edgeMap(key);
    end
    idx = unique(idx+mesh.Np);
    bdn = [bdn; idx];
end
end
